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ABSTRACT 

The magneto-rotational instability (MRI) is considered to be a promising mecha¬ 
nism to amplify the magnetic field in fast rotating protoneutron stars. In contrast to 
accretion discs, radial buoyancy driven by entropy and lepton fraction gradients is 
expected to have a dynamical role as important as rotation and shear. We investigate 
the poorly known impact of buoyancy on the non-linear phase of the MRI, by means 
of three dimensional numerical simulations of a local model in the equatorial plane 
of a protoneutron star. The use of the Boussinesq approximation allows us to utilise 
a shearing box model with clean shearing periodic boundary conditions, while taking 
into account the buoyancy driven by radial entropy and composition gradients. We 
find significantly stronger turbulence and magnetic fields in buoyantly unstable flows. 
On the other hand, buoyancy has only a limited impact on the strength of turbulence 
and magnetic field amplification for buoyantly stable flows in the presence of a realistic 
thermal diffusion. The properties of the turbulence are, however, significantly affected 
in the latter case. In particular, the toroidal components of the magnetic field and 
of the velocity become even more dominant with respect to the poloidal ones. Fur¬ 
thermore, we observed in the regime of stable buoyancy the formation of long lived 
coherent structures such as channel flows and zonal flows. Overall, our results sup¬ 
port the ability of the MRI to amplify the magnetic field significantly even in stably 
stratified regions of protoneutron stars. 

Key words: stars: magnetars - stars: neutron - supernovae: general - MHD - in¬ 
stabilities - magnetic fields 


1 INTRODUCTION 

The impact of rotation and magnetic field on the dynamics 
of core collapse supernovae is still uncertain and represents 
a promising avenue of research. In the case of fast rotat¬ 
ing progenitors giving birth to a protoneutron star (PNS) 
rotating with a period of a few milliseconds, the rotational 
energy represents an important energy reservoir that can be 
efficiently tapped if strong magnetic fields are present. In 
these conditions, the magnetorotational instability (MRI) 
has been proposed as a mechanism to efficiently amplify the 
magnetic field (A kiyama et al.|20~0 3). Its development could 
impact the explosion by converting shear energy into ther¬ 


mal energy (Thompson et al. 2005) or by generating a strong 


large-scale magnetic field which can extract the rotational 


energy and launch a magnetorotational explosion (LeBlanc 


& Wilson 

1970; 

Bisnovatyi-Kogan et al. 1976| Muller & 

Hillebranc 

lt| 19791 

Symbalisty |1984| |Moiseenko et al.| 2006 

Shibata et al. 2006; Burrows et al. 2007j; Dessart et al. 2008 

Takiwaki et al. 2009; |Takiwaki & Kotake 2011|; Winteler| 


et al.|[2012 Mosta et al.||2014 ). The magnetic field amplifi¬ 
cation by the MRI to a strength of 10 15 G could furthermore 
be a scenario to explain the formation of the most strongly 
magnetized neutron stars known as magnetars (Woods & 
Thompson 2006, and references therein). The birth of such 
magnetars with rotation periods in the millisecond range is 
a potential central engine for long gamma-ray bursts (e.g. 
Duncan &; Thompson|1 992; Metzger et al. 2011) associated 
with supernovae that have extreme kinetic energies (also 
called ”hypernovae” or type Ic BL, e.g. IlD rout et al.| |2011). 
The delayed energy injection due to the spin down of a fast 
rotating, highly magnetized neutron star has also been in¬ 
voked as an explanation for some superluminous supernovae 
like SN 2008 bi (Kasen & Bildsten 2010| Woosley| 2010 


Dessart et al.||2012 Nicholl et al.||2013 Inserra et al.|2013 ). 


Since the seminal work of Balbus & Hawley (1991), the 


MRI has been thoroughly investigated in the context of ac¬ 
cretion discs. The physical conditions in the PNS differ from 
those prevailing in accretion discs in several aspects. The 
main support against gravity is provided by pressure rather 
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2 Guilet & Muller 


than centrifugal forces, and as a consequence the rotation 
profile can be non-Keplerian. A strong differential rotation 
appropriate for MRI growth nonetheless results from the col¬ 
lapse of the progenitor core to a PNS (Akiya ma et al.| 2003 


Ott et al.|2006 ). The dependence of MRI turbulence on the 


steepness of the rotation profile was studied by Masada et ah 
( 2012 ). Another consequence of the high pressure is that ra¬ 


dial buoyancy forces driven by entropy and lepton fraction 
gradients have an important dynamical role. This stands in 
contrast to the case of thin accretion discs where only verti¬ 
cal buoyancy is dynamically relevant. The impact of buoy¬ 
ancy on the linear growth of the MRI has been studied by 


Balbus & Hawley (1994); Menou et al. (2004); Masada et al. 


(2006 2007). They showed that the diffusion of heat and lep¬ 


tons driven by neutrinos alleviates the otherwise stabilising 
effect of a stable stratification, such that the MRI can still 
grow sufficiently fast. Neutrino radiation also induces a high 
viscosity, which can limit the growth of the MRI if the initial 
magnetic field is too low (Gu ilet et al.||2015 ). 


How the non-linear phase of the MRI is affected by the 
radial buoyancy and neutrino radiation is much less known. 
The impact of buoyancy on the MRI has been studied in nu¬ 


merical simulations of ” semi-global” models (Obergaulinger 
et al. 2009; Masa da"et al.| [2015), but these simulations did 


not take into account the impact of neutrinos as a source of 
thermal diffusion and viscosity (although it has an impor¬ 
tant impact, see Guilet et al.|2015| . The pioneering simula¬ 
tions of Obergaulinger et al. ( |2009| considered a small part 
of a PNS but were not purely local in that they included 
global gradients, in particular of entropy and density. They 
investigated the impact of buoyancy driven by a radial en¬ 
tropy gradient, but artefacts induced by the radial bound¬ 
ary conditions prevented an adequate study of the impact of 
buoyancy on the non-linear phase of the MRI. Indeed, soon 
after the MRI reached non-linear amplitudes they observed 
a flattening of the entropy gradient with a strong gradient 
only present in a thin layer near the radial boundaries. The 


radially global (but vertically local) model of Mas ada et al.| 
(2015) may not suffer from such boundary artefacts (and 


could in principle describe meaningfully the evolution of the 
global radial profiles), but the radial and azimuthal reso¬ 
lution that they could afford was in our opinion much too 
coarse for an accurate study of the MRI. 

In order to avoid these shortcomings and to reduce the 
computational cost (thus allowing an exploration of param¬ 
eter space with 3D simulations), we used a more simplified 
setup. We performed shearing box simulations in the Boussi- 
nesq approximation. This approximation takes into account 
the buoyancy effects due to entropy/lepton fraction gradi¬ 
ents, but neglects other compressibility effects. As will be 
argued a posteriori, this assumption is very well justified 
since magnetic and kinetic energies remain small compared 
to the thermal energy. Our approach retains the effects of 
global gradients of angular frequency, entropy and lepton 
fraction through shear and buoyancy, but is purely local 
in the sense that the volume averages of these quantities 
are kept fixed and that the density gradient is neglected. It 
therefore retains less global effects than the ” semi-global” 


simulations by Obergaulinger et al. (2009). This seems nec¬ 


essary, however, to have clean shear-periodic boundary con¬ 
ditions which avoid numerical artefacts at the grid boundary. 


Choosing a box size smaller than the density scale-height en¬ 
sures that this local approach is approximately justified. 

The paper is organized as follows. In Section [2] we de¬ 
scribe the physical and numerical setup considered. The re¬ 
sults are then presented in Section [3] for the linear expo¬ 
nential growth phase of the MRI, and in Section [4] for the 
saturated non-linear phase that follows. Finally we discuss 
the validity of our assumptions in Section [5] and draw con¬ 
clusions in Sectionj 6 ] A linear analysis is also presented in 
Appendix 
Sections |3 


A] and [B] the results of which have been used in 
and [4] Finally, a convergence test of the simula¬ 


tions is described in Appendix [C] 


2 NUMERICAL SETUP 
2.1 Governing equations 

The simulations performed in this article are designed to 
represent a small part of a fast rotating PNS located in the 
equatorial plane at a radius of around 20 km. We describe 
the local dynamics in the framework of a Cartesian shearing 
box (e.g. Goldreich & Lynden- Bell||1965 ). The coordinates 
x , p, and z represent the radial, azimuthal and vertical di¬ 
rection, respectively, and the corresponding unit vectors are 
e x , e y , and e z . The angular frequency vector points in z 
direction, i.e. S2 = Q e z , while gravity is in — x direction, i.e. 
g — —ge x . At the location of the box considered here, the 
neutrinos are in the diffusive regime ( Guilet et al.|2015 ), and 
their effects on the dynamics can appropriately be described 
by a viscosity v as well as thermal and lepton number dif- 
fusivities. For the sake of simplicity, we assume in this work 
that the thermal and lepton number diffusivities are equal 
and denoted by x, which allows us to describe the buoyancy 
associated with both entropy and lepton number gradients 
with the use of a single buoyancy variable (see below). For 
the sake of conciseness, we refer to x as the thermal diffusiv- 
ity. We also considered the effects of a magnetic diffusivity 
77 in our study. 

As discussed in the introduction, buoyancy is described 
in the framework of the Boussinesq approximation, i.e. the 
buoyancy force is the only compressibility effect taken into 
account. The Boussinesq approximation holds if both the 
flow and Alfven velocity are much less than the sound speed 
(incompressible flow), and if the density perturbations as¬ 
sociated with the buoyancy are small. The approximation 
further assumes a uniform background density po, he. one 
neglects the radial density gradient. The validity of these 
assumptions will be assessed in Section [5. 1 | 

In the shearing box approximation, the diffusive MHD 
Boussinesq equations read in a frame rotating with the cen¬ 
tre of the box with an angular frequency £2 


dtv + v • Vr — 

dtB = 
dt0 + v-V0 = 

V • v — 
V • B = 


—vn + 

po 


47rp 0 


(V x B) x B, 


(i) 


—2f2 x v + 2qQ, 2 xe x — N 2 9e x 
V x (v x B) + 77 AS, 
v x + xA6», 

0 , 

0 , 


+ i/At) 


( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 

where v and B are the flow velocity and the magnetic field, 
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respectively. The buoyancy variable 

o = _ g_5_i 

~ N 2 po 


( 6 ) 


describes the density perturbation Sp due to the joint effect 
of entropy and lepton number perturbations, and has the di¬ 
mension of a length. The gradient of the pressure perturba¬ 
tion VII is obtained from the constraint of a divergence-free 
flow field (equation [0. The quantity 


q = —dlogfl/dlogr 


( 7 ) 


measures the local shear (we assume q = 1.25 in all sim¬ 
ulations). The square of the Brunt-Vaisala frequency N is 
defined by 

N 2 = -- 
P 

where p, S, Y e , and P are the density, entropy, electron 
fraction, and pressure, respectively. A stationary solution of 
these equations is: 


Vo 

— —qSlxe y , 

( 9 ) 

Bo 

B 0 e z , 

( 10 ) 

e 0 

= 0. 

( 11 ) 


dp I d S dp I d Y e 
dS Ip,Y g dr dY e Ip ,5 dr 


We initialised the numerical simulations with this stationary 
solution, to which we imposed random velocity and mag¬ 
netic field perturbations (as described in more detail in next 
subsection). 


2.2 Physical parameters and dimensionless 
numbers 

We chose the physical parameters of the simulations to 
represent the fast rotating PNS model studied by |Guilet| 
et al. (2015) at a radius of ~ 20km with a density po = 
2 x 10 l3 ~gcm -3 , a viscosity v — 10 lo cm 2 s -1 , and a rota¬ 
tion angular frequency SI — 10 3 s _1 . The box size was set 
to (L x , L y , L z ) — (4, 4,1) km. The radial width of the box is 
somewhat smaller than the density scale height ~ 8 km, i.e. 
the neglect of the background density gradient is roughly 


justified (see Sections 5.1 and 5.3 for a discussion). The as¬ 


pect ratio of the box, elongated in the horizontal directions 
compared to the vertical direction, was chosen to ensure 
that the box fits the most unstable parasitic modes grow¬ 
ing on 


& Goodman 

2009 

Latter et al. 2009; Obergaulinger et al. 

2009||Pessah| 

to 

0 

0 

)• 


Given these values for the box size and for the viscos¬ 
ity, the dimensionless Reynolds number (characterizing the 
importance of viscosity over the largest vertical scale in the 
box) is 


R e = = 10 3 . 

V 


( 12 ) 


According to Masada et al. (2007), the thermal diffu 


sivity x is larger than the neutrino viscosity v by a factor 
10 — 100. To investigate its impact on the results we chose 
two different values for the thermal diffusivity, namely a 
’’high” one x — 10 n cm 2 s -1 , and a ten times lower one 
X = 10 lo cm 2 s -1 referred to as ’’low diffusion” in the fol¬ 
lowing. These values of the thermal diffusivity correspond 
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to a dimensionless Prandtl number (characterizing the rela¬ 
tive importance of viscosity and thermal diffusion) 


Pr = v lx — 0.1 and 1, 


(13) 


respectively. Hence, the dimensionless Peclet number (char¬ 
acterizing the relative importance of advection and thermal 
diffusion) has a value of 


Pe = 


r 2 q 

= 100 and 1000 

X 


(14) 


for the high and low diffusion case, respectively. 

The resistivity 77 is expected to be many orders of mag¬ 
nitude smaller than the viscosity p, but numerical con¬ 
straints do not allow one to perform numerical simulations 
for such values of the resistivity. In our study, we used a 
resistivity 77 sal 2.5 x 10 9 cm 2 s _1 , which corresponds to a 
magnetic Prandtl number characterizing the relative impor¬ 
tance of viscosity and resistivity of 


Pm = p /77 4. 


(15) 


This much too small value is a compromise between the 
computing time constraints and the wish to simulate flows 
in the high ma gnetic Prandtl regime appro priate to a PNS 
(Prr 


10 1 


Thompson & Duncan 


1993 


Masada et al. 


2007). Note that this regime is unusual since the magnetic 


Prandtl number is ordinarily very small for stars, most parts 
of accretion discs and liquid laboratory metals. The high 
value of the magnetic Prandtl number in a PNS is due to 
the neutrinos inducing a very large viscosity while the resis¬ 
tivity remains very small. The dependence of the results on 
the magnetic Prandtl number will be investigated in a forth¬ 
coming publication. The corresponding magnetic Reynolds 
number , which characterizes the relative importance of mag¬ 
netic advection to magnetic diffusion (resistivity), is 

Rm = — = 4 X 10 3 (16) 

V 

in our study. 

We considered a single value of the magnetic field 
strength, whose value Bo = 1.6 x 10 13 G lies within the vis¬ 
cous regime described by Guilet et al. (2010 (see left panel 
of their Fig. 10). A dimensionless number characterizing the 
strength of the magnetic field is defined as 


p = = iq4 > 


(17) 


where va = B / \Jpo is the Alfven speed associated with 
the background vertical magnetic field. The impact of vis¬ 
cosity on the linear growth of the MRI is determined by the 
viscous Elsasser number 

E ^ v ^ = f =0 - 1 - (18) 

Since E v < 1 , our simulations belong to the viscous regime, 
where the growth rate of the MRI is significantly reduced 
by neutrino viscosity compared to ideal MRI (e.g. [Guilet 


e t al.|20l5 and references therein). Similarly, the impact of 


resistivity on the linear growth of the MRI is quantified by 
the resistive Elsasser number 


E = — 

77 T]Sl 


R n 


= 0.4. 


(19) 
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Contrary to what would be expected under realistic condi¬ 
tions, the impact of resistivity on the linear MRI is therefore 
not negligible in our setup (since E v < 1). Finally, the im¬ 
pact of thermal diffusion on the linear growth of the MRI is 
quantified by the thermal Elsasser number 


( 20 ) 


V 2 a P 

E x = = 0.01 and 0.1, 

x P 

for the high and low diffusion case, respectively. This shows 
that thermal diffusion is expected to play an important role 
in the MRI linear growth, as will be detailed in Section [3] 

Furthermore, we varied the Brunt-Vaisala frequency 
(see equation |8]), which characterizes the impact of buoy¬ 
ancy, performing a set of simulations with N 2 /PI 2 E 
{ —1.5, —1, —0.5,0,1, 2,3,4,6,8,10} both for the low and 
high thermal diffusion case. 

In order to trigger the growth of the MRI, we added 
velocity and magnetic field perturbations to the stationary 
state in the following way. We imposed Fourier modes with a 
random phase and an amplitude that follows a Kolmogorov 
spectrum. The spectrum is normalized such that the kinetic 
and magnetic energy of the perturbations amount to 1% of 
the energy in the background magnetic held. The results 
in the non-linear phase are not very sensitive to the initial 
conditions, which influence the duration of the initial expo¬ 
nential growth phase (higher amplitude perturbations lead 
to a shorter growth phase). 

Note that the simulations presented in this paper, 
meant to describe the how in a fast rotating neutron star 
at a radius of 20 km, can be rescaled to other physical sit¬ 
uations sharing the same values of the dimensionless num¬ 
bers. For example, the moderately rotating PNS model stud¬ 


ied by Guilet et al. (2015) shares the same dimensionless 


numbers at a radius of 10 km with the following parame¬ 
ters: Pi — 200s -1 , v — 2 x 10 9 cm 2 s -1 , p — 10 14 gcm -3 , 
B 0 = 7 x 10 12 G, box size (L x , L y , L z ) = (4,4,1) km, 
Prandtl numbers Pr = 0.1 — 1, and magnetic Prandtl num¬ 
ber P m = 4. The Brunt-Vaisala frequencies considered then 
range from N 2 — —6 x 10 4 s -2 to 4 x 10 5 s -2 in that situa¬ 
tion. 


2.3 Units 

Throughout the paper, we normalise our results using the 
vertical size of the domain L z — 1 km, the angular frequency 
Pi = 10 3 s -1 , and the density p 0 = 2 x 10 13 gcm -3 . Time 
is therefore measured in units of 1 ms, velocity in units of 
10 8 cms -1 , the magnetic held in units of 1.6 X 10 15 G, the 
energy density in units of 2 x 10 29 erg cm -3 , the specihc en¬ 
ergy in units of 10 16 erg/g, the energy injection rate density 
in units of 2 x 10 32 ergs -1 cm -3 , and the specihc energy 
injection rate in units of 10 19 ergs -1 g -1 . 


2.4 Numerical methods 


In order to solve the Boussinesq MHD equations @-<R 


used the pseudo-spectral code SNOOPY (Lesur & Longaretti 


2005 


2007). SNOOPY solves the 3D shearing box equations 


using a spectral Fourier method, where the shear is han¬ 
dled through the use of a shearing wave decomposition (with 
time varying radial wavevector) and a periodic remap proce¬ 
dure. Nonlinear terms are computed with a pseudo-spectral 


method using the 2/3 dealiasing rule. The time integration 
is performed using an implicit procedure for the diffusive 
terms, while other terms use an explicit 3rd order Runge 
Kutta scheme. SNOOPY is parallelized using both MPI and 
OpenMP techniques. SNOOPY has been used in the past to 


study the MRI (in the absence of buoyancy) by Lesur & 


Longaretti i 

2007), 

Longaretti & Lesur 

(2010), Rempel et al. 

(2010), and 

Lesur V Longaretti (2011 

. The Boussinesq im- 


plementation has been used to study vertical convection and 
the subcritical baroclinic instability in an accretion disc by 
Lesur & Papaloizou (2010) and Lesur V Ogilvie| (2010). 


Most simulations presented in this paper were per¬ 
formed using a standard grid resolution of (n x ,n y ,n z ) = 
(256,128,64) zones. A convergence study, presented in Ap¬ 
pendix [C] shows that this resolution is sufficient to resolve 
the dissipative scales and that the results are consistent with 
those obtained at a twice higher resolution. All our simu¬ 
lations with N 2 /Pi 2 ^ —0.5 were performed with the stan¬ 
dard resolution. At more negative values of V 2 , however, we 
observed that the more vigorous turbulence led to smaller 
dissipative scales. Therefore, we increased the resolution to 
(■ n x ,n y , n z ) — (384,192, 96) for N 2 /Pi 2 — —1 (both for P e = 
100 and P e — 1000), and to (n x ,n y ,n z ) = (512, 256,128) for 
the simulation with N 2 /Pi 2 = —1.5 and P e = 100. We note 
that in all our simulations the resolution is twice lower in 
the azimuthal direction than in the radial or vertical direc¬ 
tion. This is rather usual in numerical studies of the MRI, 
because the presence of shear makes structures more elon¬ 
gated in the azimuthal direction. We have checked that our 
results are not affected by the lower azimuthal resolution. 


2.5 Turbulent energies and energy injection rates 


We define the kinetic, magnetic, and thermal energy densi¬ 
ties as 


pf 

5‘ 

III 

to | 

"O 

o 

g 

to 

(21) 

jp — 1 d2 

mag “ 8 tt 5 

(22) 

and 


Sth = \p 0 N 2 e\ 

(23) 


respectively, where the velocity u is defined with respect to 
the equilibrium shear profile Vo according to u = v — vq. 
For brevity, we will refer to u as the velocity, and to energy 
densities as energies. 

We further denote the average of a quantity A over the 
whole computational volume by (A ). The sum of the volume 
averaged kinetic and magnetic energies obeys the following 
conservation equation 

d(-Ekin) + (£ mas ) = qQ [{Rxy) + {Mxv}] + {w) _ _ (e ^ 

(24) 

where the Reynolds stress 

B xy = pou x u y (25) 

and the Maxwell stress 

M xy = zBfB* (26) 
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are the angular momentum fluxes in the radial direction due 
to the flow and the magnetic field, respectively. The sum of 
qQ(R X y) and qQ(M xy ) gives the energy density injection rate 
into turbulent flow extracted from the shear energy of the 
flow by Reynolds and Maxwell stresses. 

W = —N 2 po9u x (27) 


is the work done by the buoyancy force per unit of time and 
per volume, i.e. (W) in equation 24 represents the energy 


exchange rate from thermal to kinetic energy through radial 
transport of entropy and lepton number caused by buoyancy. 

= pov \V x u \ 2 


24 


The remaining quantities in equation 

and e y = rj\ V x B | 2 /47r, are the energy dissipation rates due 
to the viscosity and the resistivity, respectively. 

In a time average sense, the total energy injection rate 
into turbulence, qQ [( R xy ) + (M xy )] + (IF), is equal to the 
sum of the viscous and resistive energy dissipation rates (e u + 
e v)- 


3 LINEAR GROWTH PHASE 

Let us start this section by discussing stability conditions in 
the presence of buoyancy and differential rotation both in 
hydrodynamic and MHD flows. We restrict this discussion 
to the special case considered in this paper: a local box in the 
equatorial plane of the PNS, where the angular frequency, 
entropy and composition depend only on radius. We further 
restrict the discussion to axisymmetric perturbation^] and 
to flows which are stable against the hydrodynamic shear 
instability, i.e. complying with the Rayleigh stability crite¬ 
rion k 2 > 0, where hi is the epicyclic frequency defined by 


1 d(r 4 0 2 ) 
r 3 dr 


(28) 


We first consider stability conditions in the absence of 
any diffusive process. If 


hi 2 + N 2 < 0 


(29) 


the flow is convectively unstable because the Solberg- 
Hpiland stability criteria is violated (red region in Fig. [TJ. 
In the presence of magnetic field the condition reads 

k 2 + N 2 < 4Q 2 . (30) 


This criterion is less restrictive than the hydrodynamic one 
due to the appearance of the MRI for 0 < hi 2 + N 2 < 4Q 2 
(dark blue region, and red and orange regions with blue 
stripes in the right panel of Fig. [I]). The classical MRI (i.e. 
in the absence of buoyancy effect) is obtained for N 2 = 0 
and 0 < k 2 < 4£4 2 . 

When thermal diffusion is present (but no resistivity 
and viscosity), the instability region is more extended for 
two reasons. First, thermal diffusion alleviates the stabilis¬ 


ing effect of buoyancy on the MRI (Acheson 1978 Menou 
et al.|2004 Masada et al.|2007 ). As a consequence the MRI 


1 Due to the presence of shear, non-axisymmetric perturbations 
cannot lead to exponential growth but can nonetheless be tran¬ 
siently growing. While this may not be considered as leading to 
genuine linear instability, it can be important for the non-linear 
turbulent dynamics. 


instability criterion becomes 0 < k 2 < 4£4 2 , i.e. the MRI 
unstable region now comprises a region that is stable ac¬ 
cording to equation ( |30| ) (turquoise region in the right panel 
of Fig. 0. 

Another consequence of thermal diffusion is the appear¬ 


ance of a new hydrodynamic instability (Klahr & Hubbard 
2014 Lyra|2014 ), which we refer to as oscillatory convection. 


It occurs for N 2 < 0 in a parameter regime that is stable 
according to equation (29) (orange region in Fig. [TJ . Oscil¬ 
latory convection can be considered to be an analog of semi¬ 
convection (in which thermal diffusion causes the flow to be 
unstable in the presence of an otherwise stabilising compo¬ 
sition gradient), rotation assuming the role of the stabilising 
composition gradient in the case of semi-convection. 

The parameter regime explored with our simulations is 
shown with plus signs in Fig. [l] We varied the Brunt-Vaisala 
frequency in the range —1.5 < N 2 /Q 2 <10 and fixed the 
epicyclic frequency to a value k 2 /Q 2 s= 1.5. Hence, we ex¬ 
plored three regimes: non-diffusive MRI for 0 < iV 2 < 2.5, 
diffusive MRI (i.e. allowed by thermal diffusion) for N 2 > 
2.5, and a mixed regime where MRI and oscillatory convec¬ 
tion can co-exist for —1.5 < N 2 < 0. Note that in none 
of our simulations the flow is unstable to axisymmetric 
classical (non-oscillatory) convection. 

Figure [2] shows the time evolution of the Maxwell stress 
averaged over the computational box. After a linear phase 
of exponential growth, the Maxwell stress saturates in the 
non-linear phase. In this section we analyse only the linear 
phase and postpone the analysis of the non-linear phase to 
next section. All but two simulations show a non-oscillatory 
growth, which is due to the MRI modified by buoyancy. For 
low thermal diffusion (P e = 1000; right panel), the growth 
rate varies by orders of magnitudes and becomes extremely 
small for large positive values of N 2 , up to the point that no 
growth occurs at all for N 2 /Q 2 = 10. By contrast, for high 
thermal diffusion (P e = 100; left panel), the effect of buoy¬ 
ancy on the growth of the MRI is much less pronounced. 
The MRI can still grow in a short time even for large values 
of N 2 . A good agreement is obtained between the growth 
rates obtained from our numerical simulations and the pre¬ 
dictions from the linear analysis presented in Appendix [A] 
(Fig.|3f. 

Our results are in qualitative agreement with those 
of |Masada et ah (2007 ), who showed that the MRI can 
grow unimpeded by stable buoyancy if the thermal diffu- 
sivity is sufficiently large compared to the viscosity, i.e. if 
X > i/N/Q. This condition is indeed verified for high ther¬ 
mal diffusion (P e = 100 and \l v = 10), in which case 
the growth of the MRI is only mildly affected by buoy¬ 
ancy. By contrast, the growth of the MRI is drastically sup¬ 
pressed in the case of low thermal diffusion and large N 2 /Q 2 
(P e = 1000 and x — v), i.e. if the above condition is not ful¬ 
filled. Note that, even for a high diffusivity, the MRI growth 
rates are significantly smaller than those of the ideal MRI 
(cr/Q = q/2 = 0.625). This stands in contrast to the situ¬ 
ation considered by Masada et al. (2007) and is due to the 
fact that we considered a magnetic field strength that lies 
in the viscous regime, where the MRI growth is significantly 
slowed down by neutrino viscosity (Guilet et al. 2015). 

One of our simulations with unstable buoyancy and high 
thermal diffusion (N 2 /Q 2 — —1.5 and P e = 100) showed an 
oscillatory exponential growth due to the oscillatory con- 
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Figure 1. Stability diagram as a function of the square of the Brunt-Vaisala frequency N 2 (characterizing the influence of buoyancy) 
and the square of the epicyclic frequency k, 2 (characterizing the effect of rotation and shear). The left panel applies to hydrodynamic 
flows and shows the regions unstable against classical convection (red) and oscillatory convection enabled by thermal diffusion (orange). 
The right panel applies to MHD flows and shows the region unstable against the MRI (blue) including the also convectively unstable 
regions (red and orange regions with blue stripes) and the diffusive MRI enabled by thermal diffusion (turquoise). The plus signs give 
the parameters of our simulations (those for N 2 /Q 2 = 8 and 10 are not shown). 



Figure 2. Time evolution of the box averaged Maxwell stress for simulations with different stratifications, i.e. different values of N 2 /Q 2 , 
for high (P e = 100; left panel) and low (P e = 1000; right panel) thermal diffusion. Note that the time interval is different in the two 
panels. 


vection discussed above. The transition of the fastest grow¬ 
ing mode from MRI to oscillatory convection is clearly seen 
in Fig. [3] for high thermal diffusion, where it happens at 
N 2 ~ — Q 2 in agreement with the linear prediction. In 
the regime of oscillatory convection, the growth rate pre¬ 
dicted by the hydrodynamic linear analysis presented in Ap¬ 
pendix [5] matches that of the MHD linear analysis and the 
simulations, clearly demonstrating the hydrodynamic nature 
of the instability. 


4 NON-LINEAR PHASE 


When the exponentially growing axisymmetric (i.e. y- 
independent) MRI channel modes reach non-linear ampli¬ 
tudes, non-axisymmetric secondary (or ”parasitic”) insta¬ 


bilities start to grow upon the primary modes (Goodman 
|&; Xu|[T 994| |Pessah &; Goodman||2009| |Latter et al.||2009[ 


O bergaulinger et al.||2009 Pessah 2010). This leads to the 


disruption of the channel modes, thereby ending the expo¬ 
nential growth phase. From that point on the dynamics con¬ 
sists, in most cases, of 3D statistically steady MHD turbu¬ 
lence. A somewhat different dynamics takes place, however, 
in the case of low thermal diffusion (P e = 1000) and a suf¬ 
ficiently stabilising buoyancy N 2 /Q? > 3. In that case, we 
observed the recurrent emergence of exponentially growing 
channel flows and their subsequent disruption by parasitic 
instabilities. The resulting dynamics is fundamentally non¬ 
steady due to this quasi-periodic emergence of the large scale 
channel flows. For high thermal diffusion we did not observe 
recurrent channel flows, but nonetheless noticed the devel¬ 
opment of long-lived coherent flow structures in the regime 
where the flow is stable against buoyancy. Therefore, we 
split the following discussion of the properties of the non¬ 
linear phase into two subsections. The first one we devote to 
the coherent flows (either recurrent or statistically steady), 
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Figure 3. Comparison of the growth rates as a function of the strength of buoyancy as measured by N 2 /Q 2 obtained from our numerical 
simulations (diamonds) and our linear analysis (solid lines; see Appendix |A| and |b| for high (P e = 100; left panel) and low (P e = 1000, 
right panel) thermal diffusion. The dotted line gives the prediction of the linear growth rate in a hydrodynamic flow (Appendix [bJ. Blue 
colour is used when the growth rate is negative, i.e. when it gives the damping rate. We normalized both the growth rate a and the 
Brunt-Vaisala frequency N by the angular frequency Q. 





channel 


zonal 


flow 


Figure 4. Schematic representation of the shearing box geom¬ 
etry (upper panel) and of the geometry of the coherent flows 
(lower panels): channel flows (left-hand lower panel) and zonal 
flows (right-hand lower panel). 


varying but horizontally uniform structure (i.e. flow quan¬ 
tities only depend on coordinate z), and zonal flows with 
a radially varying but azimuthally and vertically uniform 
structure (i.e. flow quantities only depend on coordinate x). 

Channel flows are related to the most unstable lin¬ 
ear modes of the MRI in case of a vertical background 
magnetic field. These modes, which have a purely vertical 
wavenumber, i.e. k — (0, 0, k z ), are known to be non-linear 
growing solutions of the MHD equations (Goodman & Xu 


1994 ). They are unstable to parasitic instabilities of Kelvin- 
Helmholtz or resistive te aring mode type (|Goodman &; Xu 


1994 Latter et al.||2009| |Pessah &; Goodman||2009| |Pessah 


2010 Obergaulinger et al. 2009), the development of which 
can stop the growth of the channel mode and trigger turbu¬ 
lence. The typical structure of MRI channel modes consists 
of a sinusoidal profile that holds for the horizontal (i.e. x 
and y) components of both velocity and magnetic field [^] 
and in the presence of buoyancy for the buoyancy variable 
0 (equation [6| , too. As can be deduced from Eqs. ( |A16| )- 
( [aT8] ), in MRI channel modes the radial (x) and azimuthal 
(y) velocity components are typically in phase |^] while the 
radial and azimuthal magnetic field components are typi¬ 
cally in anti-phase with each other ^ They are shifted also 
by a quarter of a wavelength with respect to the horizontal 
(x,y) velocity components, whereas the buoyancy variable 0 
is in phase with these components (equation |A19 ). 

Zonal flows have been observed in shearing box simula- 


while we discuss in the second subsection the time and vol¬ 
ume averaged properties of the flow during the non-linear 
phase. 


4.1 Coherent flows 

In the stable buoyancy regime ( N 2 > 0), we observe the 
appearance of either long lived (for P e = 100) or recurrent 
(for P e = 1000) coherent structures in the non-linear state of 
the MRI. Two different types of such coherent structures can 
be distinguished (Figure [4]): channel flows with a vertically 


2 In a shearing box with shear-periodic boundary conditions, the 
horizontally averaged vertical (z) component of both velocity and 
magnetic field is conserved and uniform due to the divergence 
free condition (Eqs. [4] and [5j , and has a value of zero and Bo 
respectively. 

3 In the regime where the flow is unstable to buoyancy, they can 
be also in anti-phase, i.e. the associated Reynolds stress becomes 
negative (see section 4.2 and Appendix [a|. 

4 In the regime where the flow is unstable to buoyancy, they 
can be in principle in phase too, i.e. the Maxwell stress would be 
negative. This was, however, not observed for the flow parameters 
considered in this study. 
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Figure 5. Fraction of magnetic (top row), kinetic (middle row), and thermal (bottom row) energy densities contained in axisymmetric 
structures (black + symbols), channel flows (red X symbols), and zonal flows (green circles) for simulations performed with high thermal 
diffusion (P e = 100; left column) and low thermal diffusion (P e = 1000; right column). 


tions of MRI turbulent stratified accretion discs (Johansen 


et al. 2009 Simon et al. 2012 Simon & Armitage 2014). 


They mostly possess a radially varying azimuthal velocity, 
which is associated with a radial variation of the pressure. 
In the presence of a mean vertical magnetic field, a radial 
variation of the vertical magnetic field has been reported 
recently too ( |Bai||2014| |Bai &; Stone||2014 ). 


Figure [5] illustrates the occurrence of the different types 
of coherent flows as a function of the Brunt-Vaisala fre¬ 
quency by showing the fractions of kinetic, magnetic, and 
thermal energy densities contained in axisymmetric struc¬ 
tures, channel flows, and zonal flows, respectively. To obtain 
the energy densities contained in axisymmetric structures 
we first averaged the magnetic field, velocity, and buoyancy 
variable in azimuthal (y) direction. Then we used these av¬ 
erages to compute the respective energy densities averaging 
over the box and over time. Instead, we first averaged the 
flow variables horizontally (i.e. both in azimuthal (y) and 
radial (x) direction) to compute the energy densities con¬ 
tained in channel flows. Finally, for the zonal flows, we first 
averaged the flow variables in vertical (z) and azimuthal (y) 
direction. 


Figure [5] demonstrates that the flow becomes more ax¬ 
isymmetric for increasing iV 2 , i.e. for flows that are more 
stable to buoyancy. This holds for the magnetic, kinetic, 
and thermal energies at both low and high thermal diffusion 
but is more pronounced at low thermal diffusion, in quali¬ 
tative agreement with the idea that high thermal diffusion 
somewhat alleviates the influence of a stabilising buoyancy. 
The prominence of channel flows also increases monoton- 
ically with N 2 for all three energies and at both thermal 
diffusivities. This behaviour is again more pronounced at 
low thermal diffusion, where the channel flows contain up to 
about 50% of the energy. As already mentioned the chan¬ 
nel flows have a different time behaviour at P e =■ 100 and 
P e =2 1000. At high thermal diffusion they are approximately 
steady (see Section 4.1.1), while at low thermal diffusion 
they show recurrent phases of exponential growth and dis¬ 
ruption (see Section 4.1.3). 


Zonal flows stick out most clearly regarding the kinetic 
energy (remarkably, their fraction can reach up to 90%), 
while they amount to a lesser fraction of the thermal energy, 
and an almost negligible fraction of the magnetic energy. 
Their dependence on N 2 is more complex than that of the 
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channel flows. While their importance first increases with N 2 
(up to N 2 = 4 — 6 at P e = 100, and N 2 = 2 at P e = 1000), 
it starts to decrease at still higher N 2 . As for the channel 
flows, zonal flows have a different time behaviour at high 
and low diffusion. They are approximately steady and long- 
lived at P e = 100 (Section 4.1.2), but recurrent at P e m 1000 
(Section 4.1.3). 


4-1.1 Quasi-stationary channel flows 

Figure [6] shows snapshots of the azimuthal magnetic field, 
the azimuthal velocity, and the buoyancy variable for two 
simulations at high thermal diffusion for a flow unstable 
(TV 2 /Cl 2 = — 1; left column) and stable (N 2 /Cl 2 — 10; right 
column) to buoyancy. In the latter case, the flow is much 
more axisymmetric, and we observe a channel flow struc¬ 
ture superposed to a disordered flow. 

The channel flow structure and time-evolution can be 
made more apparent by calculating horizontal averages of 
the flow variables. We show space-time diagrams of these 
averages in Fig. [t] for the simulation with N 2 /Cl 2 — 10 and 
P e — 100. A first transient phase of strong channel flow oc¬ 
curs between t ~ 100 ms and 200 ms corresponding to the 
MRI channels of the exponential growth phase described in 
Section [3] After the disruption of these channels by para¬ 
sitic instabilities and the onset of turbulence, one observes 
no clear channel flow until t ~ 300 ms. At that moment, 
another channel flow appears whose amplitude and phase is 
remarkably steady until the end of the simulation, except 
for a phase shift that takes place around t ~ 750 ms. To our 
knowledge, this is the first time that such a stationary chan¬ 
nel flow has been seen to appear in MRI turbulence (only 
recurrent channel flows were sometimes observed as will be 


discussed in Section 4.1.3). 


Let us now compare the structure of this stationary 
channel flow to that of growing linear MRI modes. From 
the space-time diagrams (Fig. [7]) it is already clear that the 
wavelength of the stationary flow (one wavelength fits into 
the box) is longer than that of the exponentially growing 
mode (a mixture of modes with n — 2 and n = 3 wave¬ 
lengths fit into the box), the latter one corresponding to the 
linear prediction Q Obviously, the quasi-stationary chan¬ 
nel flow does not correspond to the most unstable MRI 
mode. The time-averaged profiles (Fig. [7] bottom panel) 
show nonetheless interesting qualitative similarities with lin¬ 
ear MRI modes. The profiles resemble quite closely that of 
a sinusoidal function whose wavelength is equal to the ver¬ 
tical size of the box, except for the profiles of the radial 
velocity and the magnetic field which both show some in¬ 
dication of smaller scale structure on top of the dominant 
n — 1 sinusoid. Interestingly, the radial velocity, azimuthal 
velocity, and buoyancy variable are in phase, while the radial 
and azimuthal magnetic fields are in anti-phase and shifted 
by a quarter of a wavelength with respect to the velocity 
profiles. This result is in accordance with the structure of 


5 The linear analysis predicts that the n — 3 mode is the fastest 
growing mode (a/Cl = 0.066), while the n = 2 mode is the sec¬ 
ond fastest growing one with a growth rate only slightly smaller 
(cr /Cl = 0.061). Two further modes (n = 1, 4) are unstable with a 
roughly twice smaller growth rate. 



time (ms) 


time averaged profiles 



z 


Figure 7. Quasi-stationary channel flow for TV 2 /Cl 2 = 10 and 
high thermal diffusion P e = 100. The upper five panels show 
space time diagrams (from top to bottom) of the horizontally 
averaged radial velocity, azimuthal velocity, radial magnetic field, 
azimuthal magnetic field, and buoyancy variable 0. The bottom 
panel gives the vertical profiles of these quantities averaged over 
the time interval 400 ms ^ t ^ 700 ms. 
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Figure 6. Three-dimensional snapshots at time t = 553 ms of two simulations with P e = 100 for a buoyantly unstable flow with 
N 2 /Q 2 = — 1 (left column) and buoyantly stable flow with N 2 /Q 2 = 10 (right column). The three rows show the azimuthal magnetic 
field (top), azimuthal velocity (middle), and buoyancy variable 6 (bottom). Note that we adjusted the colour scales to the level of 
turbulence, i.e. they are different for the two columns. 


a linear MRI channel mode, and it suggests a link between 
the stationary channel flow and linear MRI modes. 


In view of the qualitative similarity but quantitative 
difference between the stationary channel flow and the linear 
MRI modes, we interpret the stationary channel flow as an 
MRI channel mode significantly modified by the concurrent 
turbulence, which presumably brings it to marginal stability. 
Further studies of the physics of the stationary channel flow 
and its interaction with turbulence may shed light on the 
mechanism determining the level of MHD turbulence. Using 
enlarged diffusion coefficients to roughly describe the impact 
of turbulence on the channel flow may provide the basis 
for modelling MRI saturation and could allow for a deeper 
interpretation of the numerical simulations. One could, for 
example, speculate that the level of turbulence adjusts so as 
to bring the largest scale MRI mode to marginal stability 
( Ogilvie &; Livio|2001| Guilet &; Qgilvie||2012 2013). 


4-1.2 Quasi-stationary zonal flows 

Zonal flows can conveniently be visualised with a space-time 
diagram showing the radial profile of azimuthally and verti¬ 
cally averaged quantities. Fig. [8] shows such a diagram for a 
simulation with N 2 /Pi 2 = 4 and P e — 100, where a very sta¬ 
ble quasi-st at ionary zonal flow is clearly visible. It shows 
that azimuthal velocity has by far the largest amplitude 
amounting to ~ 90% of the kinetic energy, but a modulation 
with a wavelength equal to the radial size of the box is also 
recognisable in the vertical magnetic field and the buoyancy 
variable. The modulation of the vertical magnetic field is 
shifted by a quarter of a wavelength with respect to that of 
azimuthal velocity (i.e. the maximum of the magnetic field 
corresponds to a zero of the velocity), in agreement with the 
results of Bai &; Stone| ( |2014| ). We note that the magnetic 
field zonal modulation represents a negligible fraction of the 
total magnetic energy (see Fig. [ 5 ]), but it amounts to a sig¬ 
nificant fraction of the mean vertical magnetic field. Thus, 
it may play a significant role in the flow dynamics. 

The zonal flow shown in Fig.[8]is representative of flows 
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time (ms) 


Figure 8. Space-time diagrams for a zonal flow with N 2 /PL 2 = 4 
and P e = 100. The three panels show the azimuthal velocity (top), 
the vertical magnetic field (middle), and the buoyancy variable 
6 (bottom). These quantities are averaged in azimuthal (y) and 
vertical (z) direction, and plotted as a function of time and radius 
(x). 
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Figure 9. Space-time diagrams of the azimuthal velocity for 
N 2 /PL 2 = 4 and P e = 1000. The upper panel shows horizon¬ 
tally averaged vertical profiles, where recurrent channel flows are 
clearly visible. The lower panel shows vertically and azimuthally 
averaged radial profiles, where transient zonal flows can be recog¬ 
nised just after the disruption of the channel flows. 


for Brunt-Vaisala frequencies 3 < N 2 /PL 2 < 6. In this pa¬ 
rameter regime, the zonal flow dominates the kinetic energy, 
and it is very stable and quasi-steady. At lower and higher 
values of N 2 , zonal flows tend to be more noisy and less sta¬ 
ble in that their phase changes somewhat chaotically with 
time on timescales of ten or a few tens of orbits. 



Figure 11. Three snapshots of the azimuthal velocity illustrating 
the time evolution of recurrent channel and zonal flows for a sim¬ 
ulation with N 2 /PL 2 = 4 and P e = 1000: exponential growth of a 
channel mode (upper panel), its disruption (middle panel), and 
the formation of a zonal flow after the disruption of the channel 
mode (bottom panel). The snapshots are separated in time by 
two orbital periods, and they correspond (from top to bottom) to 
times t = 1615, 1627, and 1640 ms, respectively. 


4-1.3 Recurrent channel and zonal flows 

At low thermal diffusion and in the stable buoyancy regime, 
channel and zonal flows are present too. However, contrary 
to the high diffusion case, they are not quasi-stationary. This 
is illustrated in Fig. [9] with space-time diagrams of the az¬ 
imuthal velocity for N 2 /PL 2 = 4 and P e = 1000. The ver¬ 
tical profile of the horizontally averaged azimuthal velocity, 
( v y )x,y(z ) (top panel), shows recurrent phases of channel 
mode growth and sudden disruption. The radial profile of 
the vertically and azimuthally averaged azimuthal velocity, 
( Vy)y, z (x ) (middle panel) shows that the disruption of the 
channel flow is accompanied by the formation of a strong 
zonal flow. This zonal flow then decays and another phase of 
channel flow growth begins. These successive phases of chan¬ 
nel flow growth, disruption, and zonal flow formation can 
also be seen in the left panel of Fig. |10| which demonstrates 
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Figure 10. Time evolution of various box averaged quantities showing the recurrent growth of channel and zonal flows for a simulation 
with N 2 /Pl 2 = 4 and P e = 1000. Left panel: time evolution of the kinetic energy (blue), the kinetic energy of the channel flow (green), 
and the kinetic energy of the zonal flow (red). Right panel: time evolution of the rms averages of the radial velocity (blue), vertical 
velocity (green), and the non-axisymmetric component of the vertical velocity (red). 


that the formation of the zonal flow exactly coincides with 
the disruption of the channel flow. The right panel shows 
that during the phase of exponential growth of the channel 
flow, the radial velocity (as well as the even higher azimuthal 
velocity (not shown)) is much larger than the vertical veloc¬ 
ity, consistent with a channel mode dominating the flow. 
Furthermore, the flow is mostly axisymmetric, as can be 
deduced from the fact that the non-axisymmetric part of 
the vertical velocity is significantly smaller than the total 
vertical velocity. The growth of non-axisymmetric parasitic 
instabilities disrupting the channel mode is visible in the 
sudden rise of the vertical velocity, and in particular of its 
non-axisymmetric part. The disruption of the channel flow 
is then followed by a phase of decaying non-axisymmetric 
turbulence, during which the radial velocity, the vertical ve¬ 
locity and the non-axisymmetric part of the vertical velocity 
have similar magnitudes. 


Recurrent channel flow growth and disruption by par¬ 
asitic instabilities has already been observed in numerical 
simulations of the MRI in the absence of buoyancy (e.g. 
Lesur & Longaretti 2007). Lesur & Longarettii (2007) ar¬ 


gued that this behaviour is characteristic of the MRI grow¬ 
ing close to its threshold of marginal stability, be it because 
the magnetic field is very strong, or the resistivity and vis¬ 
cosity are very high. This interpretation is consistent with 
the results of our simulations, where the MRI is brought 
close to marginal stability by the effect of buoyancy. 


Figure |TT| illustrates the phases of channel flow (top 
panel), its disruption (middle panel), and the formation of 
a zonal flow (bottom panel). The upper panel shows an al¬ 
most pure horizontally uniform channel flow, with only slight 
perturbations to the right hand side of the box. Two orbits 
later, these perturbations have grown due to parasitic in¬ 
stabilities, and the channel flow is being disrupted first on 
the right hand side of the box, and afterwards in the rest 
of the box, such that two orbits later the channel flow has 
completely disappeared and a zonal flow has formed. We at¬ 
tribute the formation of the zonal flow to the non-uniformity 
of the channel flow disruption over the radial extent of the 


box. Indeed, as the disruption of the channel flow happens 
first on the right hand side of the box (1 < x < 2), the 
Maxwell and Reynolds stresses are first quenched in that 
part of the box, while they remain high in the rest of the 
box. Before the disruption extends to the rest of the box, 
the angular momentum flux is therefore non-uniform and 
the redistribution of angular momentum creates the zonal 
flow. To our knowledge, such episodic zonal flow formation 
due to the disruption of a channel flow has never been re¬ 
ported before. Since the mechanism we propose does not 
seem directly related to the impact of buoyancy (whose role 
is mostly to bring the MRI close to marginal stability), we 
suggest that it may also take place in the classical MRI 
without buoyancy. 

The dynamics described in this section for N 2 /Pi 2 = 4 is 
characteristic of simulations with N 2 /Pi 2 > 3. In this regime, 
as N 2 increases the typical timescale between phases of 
channel flow growth becomes longer. This can be interpreted 
by the smaller growth rate of the linear channel modes (see 
Fig. §, which therefore take a longer time to grow again 
after their disruption. Note that for N 2 /Pi 2 ~ 1 — 2 and 
P e — 1000, the zonal and channel flows are more similar to 
the case of high thermal diffusion: quasi-stationary, though 
somewhat more variable than those shown in Sections 14.1.11 
and 14.1.21 


4.2 Time-averaged properties 

To characterize the properties of the non-linear state and its 
dependence on iV 2 , we performed time and box averages of 
the magnetic, kinetic and thermal energies, the energy in¬ 
jection rates, and the rms values of the three spatial compo¬ 
nents of velocity and magnetic field. For each simulation, we 
adapted the time at which the averaging is started to avoid 
the initial transient dynamics (the first exponential channel 
growth phase, sometimes followed by a transient zonal flow). 
Because of the slower growth in the stable buoyancy regime, 
some simulations had to be run for a longer time period in 
order to obtain a meaningful average over the non-linear 
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Table 1. Time averaging intervals 


Pe 

n 2 /q 2 

interval 

[msec] 

100 

^ 4 

[ 200 , 

628] 


6 

[ 400, 

1000] 


8 

[ 1000, 

2000] 


10 

[ 400, 

1000] 

1000 

< 0.5 

[ 200, 

628] 


1 

[ 300, 

1600] 


2 

[ 1000, 

3000] 


3 

[ 1000, 

3000] 


4 

[ 1000, 

3000] 


6 

[ 2000, 

6000] 


8 

[ 3000, 

10000] 


phase (see Table [I] note that 6.28 ms correspond to one or¬ 
bit). 

Figure [12] illustrates how the magnetic, kinetic, and 
thermal energies, and the energy injection rates depend 
on the Brunt-Vaisala frequency. For high thermal diffusion, 
magnetic energy and Maxwell stress decrease monotonically 
with N 2 , as one may intuitively expect. Importantly, the 
dependence gets much weaker as N 2 increases, to the ex¬ 
tent that the magnetic energy is practically independent of 
N 2 for N 2 /Q 2 ^ 3, while the Maxwell stress very slowly 
declines in this regime. As a result, the magnetic energy is 
more affected by an unstable buoyancy, than by a stable 
one. It increases by a factor of 3.3 between N 2 = 0 and 
N 2 /Q 2 = —1.5, but it decreases only by a factor 2.2 be¬ 
tween N 2 = 0 and N 2 /Q 2 m 10. In the stable buoyancy 
regime ( N 2 > 0), the thermal energy and the correspond¬ 
ing injection rate are almost independent of iV 2 , and both 
are much smaller than the corresponding magnetic quan¬ 
tities. By contrast, in the unstable buoyancy regime, they 
increase rapidly in magnitude as N 2 , up to the point where 
the thermal energy and the corresponding injection rate are 
the dominant contributions at N 2 /Q 2 = —1.5, suggesting 
that turbulence is driven then more by buoyancy than by 
the shear. Note that the increased turbulence strength in 
the unstable buoyancy regime is similar to that observed in 


the context of accretion discs by Hirose et al. (2014) and 


Hirose (2015). The geometry is nonetheless different: in ac¬ 
cretion discs the buoyancy force is mostly vertical, while in 
our setup it is radial. 

At low thermal diffusion, we found qualitatively similar 
trends for N 2 /Q 2 < 1 , but with a somewhat steeper depen¬ 
dence on N 2 , a result which agrees with the idea that a lower 
thermal diffusion allows buoyancy to act more effectively. 
The behaviour for 2 — 3 < N 2 /Q 2 < 6 is more surprising. 
All energies and energy injection rates increase in magnitude 
with N 2 , in spite of the fact that, as in the high diffusion 
case, buoyancy is actually removing energy from the flow 
(because the buoyancy injection rate is negative). This be¬ 
haviour is found in the parameter regime showing strong 
recurrent channel flows, where the energy is predominantly 
contained in axisymmetric structures (see Section 4.1). The 
increase in energy and energy injection rate is due to the 
fact that, at higher iV 2 , recurrent channel flows reach larger 
amplitudes before they are disrupted by parasitic instabili¬ 
ties. The cause of the increase of the channel flow amplitude 


at termination is unclear, but may be linked to the so far 
unknown impact of buoyancy on parasitic instabilities. In 
this respect it is interesting to note that, at high thermal 
diffusion, the amplitude at which the first phase of channel 
mode growth is terminated increases with N 2 (as can be 
seen in Fig. [5]). This may be linked to the increasing energy 
observed in the non-linear phase for low thermal diffusion. 
It also highlights the fact that the amplitude at which chan¬ 
nel mode growth is terminated is not necessarily correlated 
with the strength of MRI turbulence (which decreases with 
N 2 for high thermal diffusion). 

Figures [13] and [T4| show the dependence of the ratio of 
different flow quantities on N 2 : Reynolds stress to Maxwell 
stress (top left), kinetic energy to magnetic energy (top 
right), radial to azimuthal velocity amplitudes (bottom left), 
and radial to azimuthal magnetic field amplitudes (bottom 
right). An advantage of studying such ratios is that they al¬ 
low for a comparison of the flow properties of the non-linear 
phase to that of linear channel flows ( jPessah et al. 2006). 
Indeed, while the amplitude of a quantity (say the Reynolds 
stress) due to a channel mode is unknown if one does not 
know a priori the mode amplitude, the ratio of two quanti¬ 
ties (e.g. Reynolds stress to Maxwell stress) is well defined 
for a channel mode independently of its amplitude (see Ap¬ 
pendix 0 - Figure [l3] shows that, at low thermal diffusion, 
all these ratios vary significantly with N 2 and are in general 
agreement with the prediction from the linear modes. This 
agreement is not surprising in the regime of stable buoy¬ 
ancy where channel modes contribute a significant fraction 
of the energy, but it is surprising in the regime of unstable 
buoyancy where the fraction of energy contained in channel 
modes is insignificant. The (at least qualitative) agreement 
in the latter regime therefore suggests that linear dynamics 
leaves an imprint on the properties of turbulence even when 
channel modes do not dominate the flow. 

Interestingly, in the regime of unstable buoyancy 
(N 2 /Q 2 = — 1 and —0.5) the Reynolds stress changes sign 
and becomes negative, while it is always positive (corre¬ 
sponding to an outward transport of angular momentum) 
for the MRI in the absence of buoyancy. At P e = 1000, this 
behaviour is reproduced amazingly well by the linear analy¬ 
sis. When the Reynolds stress becomes negative, the Corio¬ 
lis force is not anymore a driver of the radial motion in the 
channel mode contrary to the classical MRI. Then the buoy¬ 
ancy force drives the radial motion instead. The magnetic 
field nevertheless still plays a crucial role in the instability 
mechanism by reducing the stabilising effect of the Coriolis 
force, because no non-oscillatory growth is possible in this 
regime in the absence of a magnetic field (see Section [3| . 
Hence, we interpret the instability in the unstable buoyancy 
regime as a mixed magneto-rotational buoyant instability, 
and the increasing ratio of kinetic to magnetic energy by 
the increasing role played by the buoyancy force compared 
to the magnetic one. 

Another interesting trend visible in Fig. concerns 
both the ratio of azimuthal to radial velocity and the ratio 
of azimuthal to radial magnetic field. These ratios increase 
tremendously with N 2 , reaching a value of 100 for the for¬ 
mer one. This may qualitatively be interpreted by the fact 
that for increasing iV 2 , the buoyancy force prevents more 
and more efficiently radial motions. In the regime of stable 
buoyancy magnetic and kinetic energy are therefore domi- 
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Figure 12. Time and volume averages of the turbulent energy densities and energy density injection rates as a function of the Brunt- 
Vaisala frequency. Panels in the left and right column correspond to high (P e = 100) and low (. P e = 1000) thermal diffusion, respectively. 
They show the magnitude of magnetic (red x), kinetic (blue +), and thermal (black squares) energies in the top panels, and energy 
injection rates due to Maxwell stress (x symbols), Reynolds stress (+ signs), buoyancy force work (circles), and sum of the three (squares) 
in the bottom panels. Positive injection rates are shown in red colour, while negative ones (i.e. energy is removed from turbulent motions) 
are shown in blue. 


nated by the azimuthal component of the magnetic field and 
the velocity, respectively. The vertical components of both 
quantities (not shown) are of the same order of magnitude 
as their radial components. 

We found qualitati vely similar trends in the case of high 
thermal diffusion^] (Fig. [l4] ) . In the regime of unstable buoy¬ 
ancy, the Reynolds stress becomes negative and the ratio of 
kinetic to magnetic energy increases. The ratio of azimuthal 
to radial velocity and the ratio of azimuthal to radial mag¬ 
netic field increase significantly with N 2 . Contrary to the 
low thermal diffusion case, these trends are not reproduced 
quantitatively, however, by the linear predictions. Still, it 
is interesting that most qualitative trends are shared be¬ 
tween the simulations and the linear modes (e.g. the ratio of 
azimuthal to radial velocity increases with N 2 ), suggesting 
that linear analysis may still provide a qualitative under¬ 
standing of these trends. 

Finally, Fig. |14| provides interesting information on the 
coherent channel and zonal flows described in Section 14.11 


6 Except for an increase of the ratio of kinetic to magnetic energy 
at high A/ -2 , which is observed for low but not for high thermal 
diffusion. 


Our simulations show that the ratios of azimuthal to ra¬ 
dial velocity and of azimuthal to radial magnetic field of 
the quasi-stationary channel flow are quantitatively differ¬ 
ent from those of the most unstable linear mode (but with 
a similar qualitative trend). This confirms our conclusion 
from Section [4. 1.1 1 that quasi-stationary channel flows, while 
related to MRI modes, have a different structure presum¬ 
ably because they are significantly modified by the con¬ 
current turbulence. The simulations further show that the 
zonal flows described in Section mu are responsible for 
the bump in kinetic energy observed in Figs. [12] and [T4| for 
3 < N 2 /Pi 2 < 6. Remarkably, at the peak of the zonal flow 
(N 2 /Pi 2 — 4 and 6) the kinetic energy is in equipartition 
with the magnetic energy. Whether this is a mere coinci¬ 
dence or a general property of zonal flows is unclear. Inter¬ 
estingly, if we subtract the contribution of the zonal flow, 
the remainin g k inetic energy (blue circles in the top right 
panel of Fig. 14) has a smooth flat dependence on N 2 , with 
no visible impact of the zonal flow (like the total magnetic 
energy). This suggests that the zonal flow does not signif¬ 
icantly impact the level of turbulence developing on top of 
the zonal flow. 
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Figure 13. Ratios of different time and volume averaged diagnostics as a function of N 2 /Q 2 for low thermal diffusion (Pe = 1000): 
Reynolds stress to Maxwell stress (top left), kinetic energy to magnetic energy (top right), azimuthal velocity to radial velocity (rms 
averages; bottom left), azimuthal magnetic field to radial magnetic field (rms averages; bottom right). Black squares give the results for 
the whole flow, blue circles when the zonal flow is subtracted, and red crosses for the channel flow. The lines show properties of the most 
unstable linear modes: discretised modes (dotted) that fit in the numerical box, and continuous modes (solid) relevant to an infinitely 
high box (which may be interpreted crudely as some weighted average of the two most unstable discretised modes). 
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Figure 14. Same as Fig. |13| but high thermal diffusion (P e = 100). 


5 DISCUSSION 


, constant shear rate p] (Section ph2| , and the local 


approxi- 


In this section, we discuss the appropriateness of some of 
the simplifying assumptions made in our study: use of the 
Boussinesq approximation (Section 5.1), the assumption of 


7 In the shearing box formalism the box integrated shear rate is 
kept constant, while its local value is allowed to vary (which is 
the case in zonal flows). 
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mat ion, in particular whether our results depend on the box 
size (Section |5.3| ). 


5.1 Validity of the Boussinesq approximation 


Because our simulations are the first attempt to apply the 
Boussinesq approximation to the problem of MRI growth 
in a PNS, it is important to assess to what extent this ap¬ 
proximation provides an appropriate description of the flow 
dynamics. Several conditions should be satisfied for the ap¬ 
proximation to be justified. Firstly, the flow velocity and 
the Alfven velocity should be much smaller than the sound 
speed, so that the flow is approximately incompressible. In 
all our simulations, this requirement holds very well, because 
- 2 /c 2 s < v\/c 2 s < 10 —4 , the sound speed c s ~5x 10 9 cms -1 


being estimated from the PNS model studied by|Guilet et al.| 
([2015). Secondly, the relative density perturbation associ¬ 
ated with buoyancy, Sp/po = — ON 2 /g , should be small. Us¬ 
ing the gravitational acceleration g ~ 6 x 10 13 cms -2 (again 
from the model of |Guilet et al.||2015| ), we find that in all 
our simulations its rms value is Sp/po < 1.5 x 10 -3 , i.e. it is 
indeed quite small. Thirdly, the Boussinesq approximation 
neglects the influence of a background density gradient, an 
assumption that is no longer justified when the size of the 
box approaches the density scaleheight. For the box size used 
in our study, the density at the radial surfaces of the box 
should actually differ from the value at the centre of the box 
by about 30%. 

The neglect of the background density gradient is by far 
the main limitation of the Boussinesq approximation in our 
simulations. Note that this limitation is common to all local 
models, even those not employing the Boussinesq approx¬ 


imation (e.g. Masada et al. 2012). It can, in principle, be 
lifted only by using a global model, or a semi-global model 


(Obergaulinger et al. 2009). So far, global models resolv¬ 
ing the MRI have only been published for two dimensional 


flows (Sawai et al. 2013 Sawai & Yamada 2014), which is 


problematic for a study of MHD turbulence. On the other 
hand, the inclusion of the density gradient in the semi-global 


models of Obergaulinger et al. (2009) renders the use of pe¬ 


riodic boundary conditions problematic, and causes impor¬ 
tant artefacts at the radial boundaries of the computational 
domain (e.g., a fast flattening of the rotation and entropy 
profiles, which prevents a meaningful study of the non-linear 
phase of the MRI). The radially global (but vertically local) 
model of [Masada et al. (2015) may not suffer from such 
boundary artefacts, but the horizontal resolution that they 
could afford is in our opinion much too coarse for an accurate 
study of the MRI. We conclude that so far no satisfying so¬ 
lution has been found to study the impact of the large-scale 
density gradient on the MRI, and that presently the use of 
local simulations in the Boussinesq approximation is proba¬ 
bly the best choice to study the impact of buoyancy on the 
MRI. 


5.2 Energy dissipation rate and timescale of shear 
profile evolution 

Another limitation of local simulations is that they cannot 
describe the evolution of the global rotation profile. In shear¬ 
ing box simulations with shear-periodic boundary conditions 


in the radial direction, the shear rate averaged over the box 
is kept constant. This is justified only as long as the rotation 
profile has not been significantly modified by the angular 
momentum transport driven by the MRI. 

One can estimate the timescale over which the rota¬ 
tion profile is significantly changed by the ratio of the total 
available shear energy to the rate at which it is extracted 
by the combined effects of Reynolds and Maxwell stresses. 
The available shear energy is estimated in the fast rotating 
PNS model studied by Guilet et al. (2015) by comparing 
the rotational energy contained in the differentially rotating 
PNS to that of a PNS in solid rotation with the same total 
angular momentum. In this way, we estimate a shear energy 
of - 4 x 10 5 ° erg. The specific rate at which energy is ex¬ 
tracted from the shear [j is between 5 x 10 16 ergs -1 g -1 in 
the regime of stable buoyancy for high thermal diffusion, and 
5 x 10 17 ergs -1 g -1 in the unstable buoyancy regime. If we 
assume that this specific energy extraction rate is uniform 
in the differentially rotating envelope of the PNS (which 
contains a mass of 1.6 x 10 33 g in the PNS model consid¬ 
ered), we estimate a total extraction rate of shear energy 
between 8 X 10 49 ergs -1 and 8 x 10 5O ergs -1 . The typical 
timescale over which the shear energy is extracted is there¬ 
fore between 500 ms in the unstable buoyancy regime and 5 s 
in the stable one, which is significantly longer than the dy¬ 
namical timescale 1/Q = 1 ms, and also in most cases longer 
than the time needed to establish the non-linear statistically 
steady state. This justifies to some extent the relevance of 
the local approach, but one should keep in mind that in some 
cases the shear profile should actually have evolved before 
the end of the simulation. 


5.3 Local approximation : dependence on box size 


One problem greatly limiting the predictive power of local 
numerical simulations, is that the result is likely to depend 
on the size of the computational domain. For the MRI in the 
absence of buoyancy and in the presence of a vertical mag¬ 
netic field, numerical simulations have indeed found that the 
level of turbulence as measured by the angular momentum 


flux is proportional to the vertical size of the box (Hawley 
et al.| |1995 ). Since the relevant box size to be adopted is a 
priori unknown, this renders very uncertain the use of lo¬ 
cal models to predict magnetic field strength and angular 
momentum transport in PNS (and casts some doubts on 
the scaling formulae given by [Masada et al. (2012), who did 
not discuss the dependence of their results on the box size). 
One may argue the density scale height to be a natural up¬ 
per limit to the scales available to MRI driven turbulence, 
such that choosing a local box of a size comparable to the 
density scale height could be an appropriate choice. Whether 
this argument actually holds can be checked only with global 
numerical simulations covering at least several density scale 
heights if not the whole PNS. While 3D simulations of such 
global models seem out of reach presently, investigating the 


8 The total rate of energy injection into turbulence can be higher 
in the regime of unstable buoyancy (up to 3 X 10 18 ergs -1 g -1 ) 
due to the buoyancy force, but this energy is not extracted from 
the shear. 
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dependence of local models on the domain size would al¬ 
ready shed more light on its potential impact on the results. 
The dependence on the box size of MRI simulations that 
include buoyancy is so far unknown, and should be explored 
in the future. 


6 CONCLUSIONS 

We have investigated the impact of buoyancy on the proper¬ 
ties of the MRI in protoneutron stars. For this purpose, we 
performed numerical simulations of a local model of a PNS 
using for the first time the Boussinesq approximation. This 
approximation has proven very useful to take into account 
the impact of entropy and lepton number gradients while 
keeping the advantages of a local model, i.e. clean shearing 
periodic boundary conditions and low computational cost al¬ 
lowing for an exploration of parameter space. We compared 
the results of our simulations to those of a linear analy¬ 
sis presented in Appendix [A] The main conclusions of our 
study can be summarized as follows : 

• The linear analysis results are confirmed by our nu¬ 
merical simulations: thermal and lepton number diffusion 
alleviates the stabilising effect of buoyancy, and allows MRI 
growth in regions of parameter space where the MRI would 
be stabilised in the absence of diffusion (Me nou et al.|2 004 


diffusion, the termination amplitude of the MRI increases 
with N 2 while the strength of turbulence decreases. This 
should be kept in mind when interpreting the results of sim¬ 


Masada et al. 2006 2007). The growth rate of the MRI is 
then largely controlled by the thermal and lepton number 
diffusion, faster diffusion leading to faster MRI growth. 

• In the non-linear phase of the MRI, stable buoyancy 
(i.e. the square of the Brunt-Vaisala frequency N 2 > 0) 
favours the appearance of large scale coherent flows: verti¬ 
cally varying (but horizontally uniform) channel flows as well 
as radially varying (but azimuthally and vertically uniform) 
zonal flows. These coherent flows can amount to a signifi¬ 
cant fraction of the total flow energies. In the case of high 
thermal diffusion, they are quasi-stationary and very stable 
over time, while they show recurrent phases of formation 
and disruption in the case of low thermal diffusion. 

• The strength of MHD turbulence (characterized by the 
magnetic energy or the rate at which energy is injected into 
turbulence) increases significantly in the regime of unsta¬ 
ble buoyancy (i.e. N 2 > 0), but it decreases only mildly if 
the flow is buoyantly stable. This general result is robust, 
i.e. it holds even when the thermal and lepton number dif¬ 
fusion is varied by a factor ten. The precise variation of 
turbulence strength with N 2 does, however, depend on the 
thermal diffusion. For high diffusion, the magnetic energy 
decreases monotonically when the flow becomes more and 
more stable against buoyancy (following the same qualita¬ 
tive trend as the growth rate). In the case of lower diffu¬ 
sion, however, the magnetic energy increases for flows that 
are very stable against buoyancy, an opposite trend to the 
growth rate which decreases dramatically. This surprising 
behaviour is observed in the regime where recurrent chan¬ 
nel flows appear, and we speculate that it might be related 
to the effect of buoyancy on parasitic instabilities growing 
on MRI channel modes. 

• Importantly, we find that the amplitude at which the 
first exponential growth of the MRI is terminated is not nec¬ 
essarily correlated with the strength of the turbulent phase 
that follows. Indeed, in the realistic case of high thermal 


ulations focusing on the MRI termination (Obergaulinger 
et al.|[2009 ), and highlights the importance of studying the 
fully turbulent phase of the MRI. 

• Other properties of MRI turbulence are strongly influ¬ 
enced by buoyancy. In particular, the ratios of azimuthal to 
radial magnetic field and of azimuthal to radial velocity in¬ 
crease dramatically in the regime of buoyantly stable flows, 
where the azimuthal components of the magnetic field and 
of the velocity therefore dominate the magnetic and kinetic 
energy. In the regime of buoyantly unstable flows, the ratio 
of kinetic to magnetic energy increases while the Reynolds 
stress changes sign and becomes negative. All of these trends 
are at least qualitatively reproduced by the properties of lin¬ 
ear unstable modes. 

The magnetic field strength estimated from the volume 
and time averaged magnetic energy in the non-linear phase 
of our simulations ranges from 3.5 x 10 14 G in the case of 
stable buoyancy with high thermal diffusion to 10 15 G in the 
unstable buoyancy regime. Such magnetic field strengths lie 
the range estimated for magnetars, and they are a few tens 
to almost a hundred times stronger than the initial mag¬ 
netic field. This provides support to the ability of the MRI 
to amplify the initial magnetic field to magnetar strength. 
Furthermore, if the MRI is at the origin of a magnetar’s 
magnetic field, our results suggest that the azimuthal mag¬ 
netic field would be significantly stronger than the poloidal 
one (in stably stratified regions of the PNS). This could have 
interesting consequences for magnetars. The observation of 
magnetar activity from neutron stars with low dipole mag¬ 
netic field ( Rea et al.|2012 2013[ 2014) indicates an internal 
magnetic field much stronger than its dipolar component. 
Furthermore, the recent detection of a phase modulation of 
the pulsed emission of two magnetars has been interpreted 
as free precession arising from a deformation of the magnetar 


due to a very strong internal toroidal field (Makishima et al 


2014 2015). The generation by the MRI of an azimuthal 


magnetic field much stronger than its poloidal component 
could be important to explain these observations. 

It has been suggested that the development of the MRI 
could impact the dynamics of core collapse supernovae ex¬ 


plosions (Akiyama et al. 2003) by converting shear energy 
into thermal energy ( Thompson et al~||2005 ) or by generat¬ 
ing large scale magnetic fields that could then drive powerful 
jet-like explosions (e.g. Moiseen ko et al.|2 006; Shibata et al. 


2006; IBurrows et al.||2007| |Dessart et al.||2008| |Takiwaki 
et al. |2009| jTakiwaki fc Kotake |2011| ). In the former case 


the effect of the MRI was parametrized by an effective tur¬ 
bulent viscosity estimate d to reach values u p to a few times 
10 13 cm 2 s -1 (see Fig. 5 of Tho mpson et al.| 2005). The effec¬ 
tive turbulent viscosity can be estimated in our simulations 
from the angular momentum flux. We found values ranging 
from ~ 5 x 10 10 cm 2 s -1 in the stable buoyancy regime with 
high thermal diffusion to ~ 5 x 10 11 cm 2 s _1 in the unsta¬ 
ble buoyancy regime. These values are two to three orders 
of magnitude smaller than those used by [Thompson et al.| 
(2005), and they would probably not provide enough vis¬ 


cous heating to affect the explosion dynamics. Concerning 
the second case, the generation of a large scale magnetic field 
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by the MRI cannot be studied in a local model such as ours, 
and has so far never been demonstrated. Though we could 
not reach definitive conclusions on this matter, our simu¬ 
lations show that stable buoyancy favours the formation of 
coherent magnetic structures on the largest scales allowed 
in our local model. The magnetic field strength we obtained 
is of the same order of magnitude as that needed to launch 
jet-driven explosions 10 15 G). In order to launch a jet, 
such a magnetic field strength is, however, needed near the 
surface of the PNS, while our simulation domain is located 
deep inside the PNS. Further numerical simulations describ¬ 
ing the outer parts of the PNS will therefore be needed to 
address this question. 


We stress that our quantitative conclusions (e.g. 
strength of the magnetic field, angular momentum flux) 
should not be considered as definitive answers due to sev¬ 
eral limitations of our study. Firstly, as discussed in Sec¬ 
tion |5.3| the results may depend on the size of the compu¬ 
tational domain (a limitation shared by all local and semi- 
global simulations). This size dependence in the presence 
of buoyancy is unknown so far and should definitely be in¬ 
vestigated in future studies. Secondly, the magnetic Prandtl 
number assumed in our simulations due to numerical con¬ 
straints (P m = 4) is much smaller than its realistic value 
(P m ~ 10 13 ). In the absence of buoyancy, MRI turbulence 
is known to be very sensitive to the magnetic Prandtl num¬ 
ber, the strength of turbulence strongly increasing with P m 


(|Lesur &; Long aretti 2007 Fromang et al.||2007 Longaretti 


& Lesur 2010). If this trend were to hold in the presence of 
buoyancy, our results should be considered as lower bounds 
on the magnetic field strength and the angular momentum 
flux. 


Many questions remain open on the MRI in the pres¬ 
ence of buoyancy. As already mentioned, it will be crucial 
to determine the dependence of our results on the box size 
and the magnetic Prandtl number. The role of the strength 
and geometry of the initial magnetic field (both of which 
are very uncertain) should also be explored. According to 
models of stellar evolution including magnetic effects ap¬ 
proximately, the initial magnetic field might be dominated 
by its azimuthal component ( Heger et al.| [2005). Is an effi¬ 
cient amplification of the magnetic field possible, if the ini¬ 
tial magnetic field is azimuthal rather than poloidal? How 
do the results depend on the strength of the initial magnetic 
field? Further numerical simulations in the framework laid 
out in this paper could answer such questions. Note also, 
that this study was restricted to the equatorial plane of the 
PNS, where the stratification is perpendicular to the rota¬ 
tion axis. It would be interesting to investigate the impact 
of buoyancy at other latitudes, thereby allowing for a range 
of angles between the radial stratification and the rotation 


Finally, this study was restricted to the inner region of 
the PNS where neutrinos are in the diffusive regime. [Guilet | 


et al. (2015) showed that in the outer parts of the PNS, 


but still below the neutrinosphere, neutrinos are in the non- 
diffusive regime on length scales at which the MRI grows. 
The influence of buoyancy on the MRI in this regime is 
unknown, and should be studied using both linear analysis 
and numerical simulations. 
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APPENDIX A: MHD LINEAR ANALYSIS 


In this appendix we give an analytical description of the 
MRI modes growing in the setup described in Section [2] 
The dispersion relation obtained in Section [AT] is equivalent 
to that obtained by Menou et al. (2004) and Mas ada et al.| 
(2007), but restricted to the particular setup considered in 


this paper, i.e. for a purely vertical magnetic field, the de¬ 
scription of buoyancy effects by a single buoyancy variable 
(relevant to the particular case where the thermal and lepton 
diffusivities are equal), and considering only the dynamics 
in the equatorial plane. We repeat the derivation of the re¬ 
spective dispersion relation below, because for our setup the 
derivation is easier to follow, and because intermediate re¬ 
sults (relations between the different velocity and magnetic 
field components of the modes) are subsequently used in 
Sections IA2I and IA31 


Al dispersion relation 


For simplicity, the analysis is restricted to MRI channel 
modes whose wave vectors are purely vertical (i.e. having 
only a z component), because these are the fastest growing 
modes in the presence of a purely vertical background mag¬ 
netic field. Perturbations from the stationary background 
flow are assumed to have a time and space dependence 
SA = Re[5Ae at + ikz ] for A € {u x , u v , u z , B x , B y , B z , 0}, 
where cr, k, and 5A are the growth rate, the vertical wavevec- 
tor, and the complex amplitude of the mode, respectively. 
For conciseness, we drop the “symbol when referring to the 
complex amplitudes in the rest of this section. The equa¬ 
tions governing the time evolution of such perturbations are 
obtained from Eqs. 0 - 0 - They read 


crSux 

= a Bz ik5B x + 2Q5u v N 2 50 
4np 0 

k 2 

aSuy — SSu x 

B z 2 

— - ikSB v — 2 Q5u x — k vbu v 

47rpo 

(A2) 

Su z 

= 0 

(A3) 

crSB x 

= B z ik5u x — k 2 rjSB x , 

(A4) 

a 5 By 

— B z ik5uy SSB X k tjSB y , 

(AS) 

5B Z 

= o, 

(A6) 

cf50 

— Su x — k 2 x50, 

(A7) 


where S = qQ is the shear rate with q given by equation 0. 
Note that these equations are valid for any amplitude of the 
perturbations, as no linearization had to be done to obtain 
them. The channel mode solutions that we obtain below 
are therefore non-linear solutions, just like the classical MRI 


channel modes in the incompressible limit (Goodman V Xu 


1994). Next, we define 


= a + k v, 


(A8) 
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with 


Gy - G 

+ k 2 7], 

(A9) 

«5 

(Z4 

= 1, 

= k 2 (x + 2v + 2rj) , 

— N 2 *f k 2 + 2 k 2 v A 

+k 4 [2 x(r] + ^) + + v 2 + g 2 ] , 

cr x = cr 

+ 

to 

>< 

] 

] 

] 

] 

] 

(A10) 

CL2 

= N 2 k 2 (y + 2g) + ft 2 /c 2 (x + 2g) 

+2 k 4 v\(y + g + x) 

and rewrite Eqs. (Al), (A2), (A4), (A5), and (A7) 

as 


-\-k 6 [x(^ 2 + g 2 + 4?7 v) + 2vg(g + v 

Gi/SlL x 

= A ik5B x + 2QSu y N 2 50, 

47T/90 

= a Bz ikSBy <2Q S )5u x , 

47 t po 

(AH) 

a\ 

= k 8 [g 2 v 2 + 2x(gv 2 + vg 2 )\ 

-\-N 2 k 4 (g 2 + 2gu) + K 2 k 4 (g 2 + 2gx) 

Gjy 6lJjy 

(A12) 


+k 2 v\ [2k 4 (vg + v\ + gx) + 

Gy 6 B X 

— B z ik6u x , 

(A13) 


+k 2 v 2 a + k 2 — 4U 2 ] , 

G y 6 -By 

— B z ik6u y S6B X , 

(A14) 

do 

7 10 2 2 . a t -2 7 6 2 | 2/6 2 

= k XO v + N k vg -\- k k XO 

g x 60 

— 6u x . 

(A15) 


+k 2 v A [2k 6 vgx + N 2 k 2 g 


Using Eqs. (All) and (A13)—(A15), one can express all 


non-vanishing variables as a function of the radial velocity 
amplitude 


+{k 2 v\ + « 2 - 4 n 2 )k 2 x]. 


(A22) 

(A23) 

(A24) 


(A25) 


(A26) 


(A27) 


In order to obtain the growth rates, this fifth order polyno¬ 
mial is solved numerically using a Laguerre method. 


1 


r * . N 2 k 2 v 2 A \ . 
° U V — 7)7T ( Ov H I ) 6u xr 

Zi\ L \ G x Gy 


5B X = B- 


ik5u x 


5 By — ( cr „ H-b 


60 = 


N 2 k 2 v\ — 2QS\ B z ikSu x 


J G y 


6u x 


(A16) 


(A17) 


(A18) 


(A19) 


Using equation ( |A12| , we then obtain the dispersion relation 

cr x + k 2 v\) 2 + K (cr 2 + k 2 v\) ~ 4Q 2 k 2 

+N 2 g v {g v g v + k 2 v A ) = 0, (A20) 

where k is the epicyclic frequency defined in equation (28}. 


The more general equation (31) of Masada et al. (2007) re¬ 
duces to our dispersion relation in the case of a purely ver¬ 
tical magnetic field and wavevector, provided the thermal 
and lepton number diffusivities are equal (in which case our 
Brunt-Vaisala frequency should be identified with the sum 
of two frequencies of Masada et al. (2007) that are related 


to the entropy and lepton number gradient, respectively) or 
either the lepton number gradient or the entropy gradient 


vanishes. Equation (A20) is also equivalent to equation (13) 
of Menou et al. (2004), when applied in the equatorial plane 


to a purely vertical magnetic field and wavevector. 

The dispersion relation is a fifth order polynomial in cr, 
which can be written in the form 


CLscr 5 + a^G 4 + a3cr 3 + a2cr 2 + mcr + ao = 0, (A21) 


A2 Reynolds and Maxwell stresses 

The vertical average of the Maxwell and Reynolds stresses 
associated with a channel mode can be computed using the 
complex amplitudes of the velocity and magnetic field per¬ 
turbations. The vertical average of the product of two quan¬ 
tities / and g is obtained from their complex amplitudes / 
and g through 

rL z /2 


(fg) = r 


U 


-L z /2 


fgdz = -R e(fg*), 


(A28) 


where g* is the complex conjugate of g. Here, we 
have used the relation Re(z\ ) Re (22) = (z\ + z*)(z 2 + 
z %)/4 = [Re(ziz$) + Re(z\Z 2 )] /2, for two complex num¬ 
bers z\ and Z 2 , and the fact that the vertical average of 
Re (^fge 2at+2,lkz ^j vanishes (for non vanishing k). 


Using equation (A16), the vertical average of the 
Reynolds stress reads 


R x 


N 

— Re ( <7 U H-b 


k 2 v\ 


po 


\6u x 


4U 


(A29) 


For an unstable mode (cr > 0) and N 2 ^ 0, the Reynolds 
stress is always positive, which is usually the case for MRI 
channel modes and turbulence (e.g. Pessah et al. 2006). 


When N 2 ^ 0 on the other hand, the buoyancy can po¬ 
tentially make the Reynolds stress take negative values. 


M x 


Using Eqs. (|A17|) and (|A18|, we obtain the vertical av- 

k 2 v 2 A -2QS\ 


erage of the Maxwell stress 

r 7 2 2 
k v A 


— —Re 


N 


H-b 
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xpo 


| Sv 


4U 


(A30) 


Using the dispersion relation, this expression can be rewrit¬ 
ten as 


(Jn 


N 2 k 2 v\ 
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(A3l) 

The same remarks as for the Reynolds stress apply for the 
Maxwell stress. However, note that the term proportional 
to ft 2 makes it less likely for the Maxwell stress to become 
negative, because it would require more negative values of 
N 2 than for the Reynolds stress. 


A3 Kinetic and magnetic energies 

In a similar way as for the Reynolds and Maxwell stresses, 
one can obtain the vertically averaged kinetic and magnetic 
energies 


Ek — Re 


f a^ + N 2 /a x +k 2 v 2 A /ar,^ 2 

4n 


po\Su. 


■X 2 , 


(A32) 


Em — Re 


1 + 


cf u + N 2 /a x + ( k 2 v 2 A — 2QS)/a v 


4Q 


rv v j 4 | _ 12 

X R^I) po]SUxl ' 


(A33) 


APPENDIX B: HYDRODYNAMICAL LINEAR 
ANALYSIS 

The dispersion relation in a hydro dynamical flow can be 
obtained from the MHD dispersion relation by setting the 


Alfven velocity to zero in equation (A20) 


a x (cr 2 + ft 2 ) + N 2 cr v — 0. 

(Bl) 

This third order polynomial in a may 
form 

be written also in the 

a3<r 3 + d2(J 2 + aicr + ao = 0, 

with 

(B2) 

a 3 - 1, 

(B3) 

d2 — k 2 x + 2 k 2 u, 

(B4) 

ai = k^v 2 +2k A vx +k 2 + N 2 , 

(B5) 

ao = k 2 x (k 4 v 2 + ft 2 ) + k 2 vN 2 . 

(B6) 


In order to obtain the growth rates shown in Fig. [3] we solve 
numerically this third order polynomial using a Laguerre 
method. 

An interesting approximate analytical solution of the 
dispersion relation can be obtained in the inviscid limit by 
performing an asymptotic expansion for N 2 <C ft 2 . The ze¬ 
roth order solution is a simple epicyclic oscillation of fre¬ 
quency ft and vanishing growth rate. Performing the expan¬ 
sion up to first order provides the complex growth rate for 
small N 2 

N 2 1 N 2 k 2 X 


1 + 


2 (k 4 x 2 + ft 2 ) 


(B7) 


2 (k 4 x 2 + ft 2 ) 

where the imaginary part of the growth rate is the oscillation 
frequency and the real part is the growth rate. The fastest 
growing mode has the following wavenumber and complex 
growth rate 


^max — \f ft /Xi 


(B8) 


1 + 


TV 2 ' 
4 Hi 2 


A 2 
4 Hi ’ 


(B9) 


We found this analytical solution to be a rather good approx¬ 
imation in the inviscid limit even for N 2 of the order of ft 2 . 
The maximum g rowth rate that we obtain is the same as the 
one derived by Lyra (2014) for N 2 <C ft 2 and a purely verti¬ 
cal wavevector (although these authors considered a thermal 
relaxation for a given thermal time rather than the diffusion 
approximation that we used here to describe length scales 
longer than the neutrino mean free path). 


APPENDIX C: CONVERGENCE TEST 

We performed a convergence test to check whether the spa¬ 
tial resolution employed in our simulations was sufficient to 
resolve the dissipation scales, and hence whether the results 
are unaffected by the limited grid resolution. The spectra of 
the energy and dissipation rates are shown in Fig. |B1| for 
N 2 — 0 and P e = 100, comparing a run performed with our 
standard resolution [n x , n y , n z \ = [256,128,64] and a run 
performed at a twice higher resolution. The energy and dis¬ 
sipation spectra are in good agreement for kL z l2ii < 20, i.e. 
in a regime where most of the energy is located and most of 
the dissipation occurs. The dissipation rate for kL z /2 ty > 20 
(where the spectra start to differ) is at least an order of mag¬ 
nitude smaller than its peak value. Thus, these unresolved 
scales have little impact on the dynamics. 

When varying the Brunt-Vaisala frequency, we observed 
for N 2 < 0 the peak of the dissipation rate shifted to smaller 
scales due to the enhanced turbulent activity. We therefore 
increased the resolution in order to ensure that the peak of 
viscous, resistive and thermal dissipation is well resolved. 
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